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.1  bstract 

T  his  paper  traces  the  development  of  the  understanding  of  contaminant  transport  in  an 
aiptifer  for  a  radially  symmetric  region.  It  presents  the  progression  of  ideas  and  equations  lead¬ 
ing  to  the  equation  set  describing  the  advective/dispersive  mechanisms  coupled  with  rate-limited 
adsorption.  This  equation  set  is  converted  to  a  numerical  scheme  in  the  Fortran  language  u-ing 
finite  elements  and  finite  differences.  The  resulting  model  is  tested  against  the  Laplace  t  rau-h  am 
solution  to  the  same  equation  set,  and  severs,  graphs  are  presented  detailing  the  rompari.-«  at 


NTMKHIC’AL  MODELING  OF  CONTAMINANT  T HANSFOH  ! 


with  rate-limited  sorti  ion/dksohp no.x 

IN  AN  AQl’II  KR 

/.  Introduction 

I  I  Tin  in (  unit  dual s 

'J  Years  nf  accident al  or  > l<  1 1 1 >> •  r.« t •  •  dumping  of  petroleum  products  au.|  hazardous  wa.-r*  iim 
lerials  has  resulted  in  polluteil  groundwater  I  he  Department  of  Defense  is  engaged  m  a  m.i."h<' 
program,  known  as  the  Installation  Restoration  Program,  to  cleanup  these  polluted  environin' nt> 
on  military  installations,  and  the  I'nvironmenta!  Protection  Agency  is  funding  research  to  help  in 
the  cleanup  efforts.  These  cleanup  efforts  often  involve  drilling  wells  in  t  he  vicinity  of  I  he  coni  am 
mated  area  and  evacuating  water  until  the  concent  rat  ion  of  oollutant  is  reasonably  low  i  d'J'.t  i. 
j  i pi.  ( :,2:t;:to). 

I  h * •  planners  at  t liese  remediation  sites  use  tools  to  est innate  how  long  the  wells  must  oj,.  rate 

base, |  on  some  initial  sample  data  and  sound  judgment  ,  These  tools  are  frequently  numerical 

*•  y 

simulations  of  the  advective/dispersive  transport  equation  that  are  fine-tuned  for  each  type  .  I 
pollutant  and  aquifer  geometry  To  assist  the  modelers,  studies  have  been  conducted  in  laU-rat.  re  - 
to  determine  the  fate  of  pollutants  in  vara  ms  soils  (.'id.  a  I  Many  recent  elTorts  at  tempt  t< .  .  •nq  in 
\ ai  ions  mathematical  models  to  laboratory  data  and  In  Id  data  (*’>.  21.  IS.  .'ill.  1  ti.  lab  ■  a  •  \ ■  n  t • 
simulated  data  (d~.  '.))  lies'  investigations  asstws  the  impact  of  adsorpt  ion  <•(  organic  pollutant' 
in  different  types  of  soil  lliis  information  is  used  to  upgrade  tin  accuracy  of  tin  models  .m  l  o> 
adjust  tin-  methods  employed  m  the  mo.h  Is  and  tl;e  remediation  efforts. 

I  he  numerical  modeling  of  contaminant  transport  In  polluted  groundwater  is  an  impoitant 
part  of  the  remedial  ton  |  he  knowledge  gained  lay  running  t  hese  models  can  be  used  by  plan  in  i~  to 
make  p  reel  ict  nuts  about  the  cost  and  dur.it  ion  of  cleanup  elforts,  and  to  gain  furt  her  under >t  an  ding 


of  tin"  physical  processes  involved  in  the  transport  of  foreign  materials  in  the  aquifer  its'  If  1  b- 
ad  vert  ion  /  dispersion  equation  has  traditionally  been  used  to  model  the  transport  of  contain. n 
It  uses  a  retardation  factor  to  account  for  the  chemical  sorbing  to  the  aquifer  materials,  imp!;,  iuj 
the  local  equilibrium  assumption  (LKA),  where  equilibration  between  sorbed  and  aqueous  pi,  i-  s 
occurs  instantly  (2). 

With  this  thesis  I  propose  a  numerical  solution  of  the  mathematical  model  of  the  standard 
transport  equation  (2)  modified  to  incorporate  diffusion  into  and  out  of  regions  with  ini'n-  ■  1  ■  i ! ■ 
’••ater.  These  regions  of  immobile  water  may  contribute  to  rate-limited  sorptiou/desorpt  i« -a  1  1 
contaminant,  and  have  a  significant  impact  on  aquifer  remediation.  Sand  and  gravel  aquifer-  .u- 
the  only  ones  considered  in  t  his  research .  and  the  unsat  u  rated  zone  is  ignored  (the  uusat  urat-  d  :  t  ■ 
is  the  zone  just  above  the  water  table).  The  mail-  -terns  of  interest  in  this  effort  arc  the  mod-  bug 
of  the  breakthrough  curve  as  the  contaminant  is  pumped  out  of  the  aquifer,  the  retardation  --I 
the  contaminant  due  to  sorption,  the  contribution  of  the  immobile  zone  diffusion  to  the  -tailing' 
of  the  breakthrough  curve,  and  tlx-  effect  of  pulse  pumping  (variation  in  the  pumping  ratei  a 
remediation. 

The  concentrations  in  the  water  extracted  from  a  well  during  pump-and-t real  remediate  -a 
is  governed,  in  part,  by  the  transfer  of  the  chemicals  from  a  reservoir  of  contaminant  (sorbed  -  r 
otherwise)  to  the  flowing  water  (22:t>32). 

The  goals  of  this  research  an-: 

1.  Create  a  numerical  model  that  can  he  run  on  a  small  computer.  'The  model  will  lie  ha.-  1  mi 
the  transport  equation  and  will  include  the  physics  of  rate-limited  sorpt  ion/desorpt  i- >u  in  a 
heterogeneous  aquifer. 

2.  i  sing  the  model,  demonstrate  the  breakthrough  curve  and  the  observed  tailing  for  > 1 1 ; t  if  ■  !•■ 
choices  of  the  input  parameters.  This  will  partially  validate  the  numerical  s« -hit  i •  ■  1 1  .  f  th- 
diffusion-limited  transport  mechanism,  and  give  insight  into  the  effect  of  rate -limited  ad- op¬ 
tion  on  the  fate  of  these  pollutants. 

•  i.  Demonstrate  how  a  variation  in  the  pumping  rate  (pulse  pumping)  might  he  used  advanta¬ 
geously  by  those  involved  in  the  remediation,  and  how  slow  desorption  will  effect  the  cmiiaiii- 
inant  concentrations  after  the  pumping  stop-. 

1  Cain  a  practical  understanding  of  finite  difference  and  finite  element  methods  use-1  in  nmd- 
eling. 


/.  2  Background 


The  simplest  contaminant  trans|iort  model  incorporates  advection,  dispersion,  and  pi 1 1 1 1 >- 
rium  sorptiifi  to  simulate  the  movement  of  a  pollutant  in  an  aquifer.  This  sin  pie  model  mn\  l.e 
suitable  when  relatively  low  levels  of  contamination  are  not  of  significance.  However,  when  very  low 
levels  of  contamination  are  significant,  such  ..es  during  remediation  of  aii  aquife'  ontaminated  with 
an  organic  pollutant,  the  simp,  advective/dispersive  niodel  is  .nsufificient  for  demonstrating  the 
tailing  that  has  been  found  in  field  and  in  laboratory  experiments  (19).  One  explanation  fo.  this 
tailing  (whereas  the  advective/dispersive  model  predicts  a  virtually  symmet  ric  breakthrough  curve) 
is  the  existence  of  areas  in  the  aquifer  where  the  water  doesn't  flow.  The  contaminants  diffuse  into 
these  immobile  regions  due  to  the  existence  fa  concentration  gradient  This  is  further  exacerbated 
by  there  typically  being  a  larger  porosity  in  these  low  permeability  immobile  zones  (22:62  1).  As 
the  well  pumps,  and  the  main  mass  of  pollutant  ;s  adverted  away,  the  gradient  changes  and  the 
chemicals  legin  t„  diffuse  out  of  these  immobile  zones  back  into  the  moving  water  (the  mobile 
region).  As  the  chemicals  diffuse  out  of  these  immobile  regi  ms,  the  resultant  liux  produces  only  a 
slow  decrease  in  the  level  of  contaminant,  extracted  from  the  well,  or  experience!  in  the  laboratory 
(the  tailing). 

Coats  and  Smith  proposed  a  mathematical  model  to  account  for  this  phenomenon.  I  hey 
postulated  the  existences  of  immobile  zones  due  to  dead  end  pores  (10).  In  their  model,  the 
contaminant  moves  between  the  mobile  and  immobile  regions  at  a  rate  proportional  to  the  gradient 
of  the  concentrations  in  the  two  regions.  Their  method  was  avoided  for  many  years  because  of 
computational  difficulties  (21:276). 

/../  Statement  of  the  Problem 

Colt?.  (14)  (11!)  presented  an  analytical  model  which  includes  Fickian  diffusion  in  immobile 
zones  It  uses  diffusive  exchange  between  the  two  regions  (mobile  zone  and  immobile  zone)  ,  aid  uses 
t  lie  advective/dispersive  equat  ion  to  model  transport  in  the  mobile  region,  t  hits  allowing  for  the  si  uv 
sorption  and  desorption  of  contaminant  between  mobile  and  immobile  zones.  This  mathematical 
model  assume*-'  axial  symmetry  of  the  contaminant  from  an  extraction  well  out  to  sonu  arbitral) 


distance.  The  initial  contaminant  concentration  from  the  well  out  is  assumed  constant,  dropping 
to  zero  at  the  border  of  the  contaminated  region.  A  model  that  allows  for  more  general  boundary 
and  initial  conditions  is  one  step  toward  a  model  that  may  be  used  by  planners  at  a  cleanup  site. 


l.J  Scope  of  the  Research 

This  research  details  the  creation  of  a  numerical  model  based  on  the  physical  principle  of 
mobile  and  immobile  zones  with  the  possibility  of  Fickian  diffusion  in  the  immobile  zone.  I  he 
numerical  model  will  be  compared  with  the  analytic  solution  of  Goltz  (13)  to  confirm  the  correct  ness 
of  the  model. 


/.  T  Assumptions 

There  are  several  assumptions  in  the  formulation  of  the  model: 


•  With  any  flow  rate  there  will  be  a  drawdown  (see  Figure  1.1).  In  the  model,  drawdown  o! 


//"/////////////"/ 

Idealized 


Actual 


FigUi  ’1.1.  Drawdown  assumption  around  a  fully  [tenetrating  well  in  an  aquifer 

the  height  of  the  aquifer  due  to  the  pumping  is  ignored. 

•  The  model  assumes  a  homogeneous  aquifer  material. 

•  The  height  of  the  aquifer  is  constant,  the  flow  is  not  fractal,  and  the  background  seepage  rate 
is  negligible  when  compared  with  the  water  movement  due  to  pumping. 
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•  The  contamination  has  a  radial  symmetry  and  is  throughout  the  vertical  extent  in  the  aquifer 
(see  Figure  1.2).  In  reality  the  concentration  distribution  would  depend  upon  how  the  aquifer 


Figure  1.2.  Radial  Symmetry  of  contaminant  profile  versus  Reality 


was  contaminated.  It  is  also  assumed  that  the  concentration  is  limited,  and  no  further  con¬ 
tamination  takes  place  i.e.,  no  external  sources  or  sinks  of  pollutant. 

•  The  molecular  diffusion  in  the  mobile  region  is  much  smaller  than  the  dispersion  caused  by 
the  velocity  induced  by  the  pump  (33). 

1.6  Materials  and  Support  Requirements 

Successful  completion  of  the  goals  of  this  thesis  required  the  availability  of  user-frieudl\ 
computer  hardware  and  a  standardized  computer  language  accepted  by  the  scientific  community. 
Both  the  Sun  workstation  and  various  high-end  personal  computers  were  considered  as  candidate 
hardware  systems.  The  Fortran  language  was  the  obvious  choice  for  this  scientific  application 
because  of  the  availability  of  optimized  compilers  and  numerous  well-tested  numerical  recipes  (27). 
The  Sun  workstation  w'as  readily  available  and  liad  the  necessary  compiler  and  graphics  utility 
( Mathematica),  but  lacked  the  IMSL  routines.  I  chose  the  Amiga  over  the  other  available  high-end 
personal  computers  because  of  its  system-level  multitasking,  reasonable  price,  speed,  and  odor 
graphics.  I  chose  the  Absoft  AC/FORTRAN  compiler  (1)  because  it  implements  the  entire  ANSI 
X-3.9  1978  standard  (commonly  known  as  Fortran-77)  (23).  As  it  happened,  the  Absoft  compiler 
compiled  the  source  code  considerably  faster  on  the  Amiga  than  the  compiler  on  the  Sun  3/r>0! 


/.  /  Overview 


The  overall  thesis  effort  consisted  of  four  key  phases.  We  begin  with  a  survey  of  recent 
literature  in  the  field  of  modeling  sorption  processes  during  contaminant  transport,  and  develop  the 
mathematical  expressions  implemented  in  the  numerical  model.  Then  we  develop  the  numerical 
schemes  used  to  transform  the  system  of  partial  differential  equations  into  a  form  suitable  for 
programming  on  a  computer,  and  briefly  examine  the  numerical  properties  of  stability,  consistency, 
and  convergence  for  the  numerical  methods.  And  finally,  we  present  the  results  of  sample  executions, 
and  compare  them  with  the  computer-based  analytical  solutions  of  Goltz  (13). 

The  main  part  of  the  thesis  is  presented  here,  with  the  source  code  available  as  an  appendix 
from  the  Mathematics  and  Computer  Science  Department  of  the  Air  Force  Institute  of  Technology. 
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II.  Mathematical  Formulation  of  the  Problem 


2.1  Introduction 


Cleanup  procedures  at  sites  where  toxic  materials  have  contaminated  the  groundwater  * <('t <-u 
include  pumping  the  contaminated  water  out  of  the  aquifer  until  the  level  of  contamination  is 
acceptable  (9:1329).  The  models  presently  in  use  to  help  in  the  aquifer  cleanup  efforts  do  not 
account  for  the  possibility  of  rate-limited  sorption/desorption  and  molecular  diffusion  in  regions 
of  low  permeability  (13,  22).  Due  to  rate-limited  desorption,  the  toxic  materials  r  ay  persist 
longer  and  at  a  much  higher  concentration  than  predicted  by  models  which  assume  equilibrium 
sorption  /  desorption . 

The  following  sections  examine  the  efforts  of  researchers  to  develop  a  mathematical  model  to 
account  for  more  of  the  known  physical  processes  in  contaminant  transport.  These  mathematical 
models  necessarily  precede  the  development  of  computer  based  solutions. 

2.2  The  Advective/ Dispersive  Model  and  Linear  Sorption 

The  traditional  way  to  model  contaminant  transport  is  with  advection,  dispersion,  and  some 
exchange  of  contaminant  with  the  solid  material  in  the  aquifer.  Equation  2.1  shows,  in  cylindrical 
coordinates,  the  mass  balance  typically  used  to  account  for  these  processes. 


d£  pdS_ 
dt  +0  8t 


l_d_ 

r  dr 


f2.1) 


where 


C(r,t)  :  contaminant  concentration  in  the  mobile  zone  [M/L3] 
r  =  radial  coordinate  [L] 
t  =  time  [T] 

S(r.t)  =  sorbed  contaminant  [unitless] 

V’(»  )  =  seepage  velocity  [L/T] 
p  =  bulk  density  of  aquifer  material  [M//.3] 

0  —  aquifer  porosity  [unitless] 
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The  term  D  =  A/|V(r)|  +  D*  is  the  hydrodynamic  dispersion  coefficient  [L2/T],  Ai  is  the 
longitudinal  dispersivity  of  the  porous  medium  [£],  (Cf?-)!  is  the  magnitude  of  the  seepage  velocity 
[L/T],  and  D*  is  the  molecular  diffusion  coefficient  [L2 /T].  The  seepage  velocity  is  given  by 


V(r) 


-Qw 

2tr  r  H  6 


(2.2) 


where  Qw  is  the  pumping  rate  at  the  well  [Z,3/T],  H  is  the  mean  height  of  the  aquifer  [L],  and  0  is 
the  porosity  or  volumetric  moisture  content  of  the  aquifer  [unitless]. 


The  molecular  diffusion  is  usually  taken  to  be  very  small  compared  to  the  dispersion,  and  the 
term  is  dropped.  This  is  certainly  true  near  the  well  when  the  well  is  pumping,  but  becomes 
less  accurate  as  we  get  further  from  the  well,  or  as  the  flow  rate  at  the  well  decreases. 

The  many  models  currently  in  existence  differ  primarily  on  how  the  —  term  is  represented 

at 

(6:34)  .  Nonlinear  sorption  isotherms  are  typically  represented  with  the  Freundlich  equation:  .S  = 
h'd  Cn  (6)[35].  The  standard  procedure  assumes  equilibrium  and  a  linearity  between  the  time 
rates  of  change  of  the  sorbed  and  aqueous  solvent,  viz.,  5  =  Kd  C,  where  I\d  is  the  distribution 
coefficient  [L3/ M]  and  n  =  1.  Substituting  this  into  Equation  2.1  produces: 


R 


dC(r,t )  1  d 


dt 


dr 


pl<< 


r  D 


dC(r ,  t ) 


dr 


-  T(r) 


dC(r,  t) 
dr 


(2.3) 


Where  R  —  1  +  is  the  retardation  factor  [unitless].  This  will  account  for  some  retarding 

V 

of  the  pollutant  as  it  temporarily  adheres  to  the  materials  in  the  aquifer.  This  is  known  as  the 
local  equilibrium  assumption  (or  LEA)  (3). 


2. -1  Sorption  Kinetics 

One  of  the  very  early  attempts  to  mathematically  model  diffusion  processes  in  adsorption 
columns  is  that  by  Lapidus  and  Amundson  (20).  In  their  1952  paper  they  propose  a  diffusion 
process  as  the  cause  of  slight  deviations  between  predictions  and  laboratory  results. 


The  asymmetry  in  the  effluent,  concentration  profiles  from  laboratory  experiments  with 
columns  of  soil  led  Coats  and  Smith  (10)  to  postulate  the  existence  of  immobile  zones  within 
the  aquifer  material.  They  proposed  a  differential  capacitance  mathematical  model  equivalent  to 


that  of  Lapidus  and  Amundson,  with  diffusion  into  small  stagnant  volumes  of  water.  Their  equal  i< >ii 
set  in  one  spatial  dimension  is  given  here  ( 10:76): 
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ac  „  as 


(i 


A)£ 

’  at 


a2c  ac 

—  D-—-: —  V  —  mobile  region 

dx-  dx 

—  K  (C  —  S)  immobile  region 


(2.4) 


Where  S  is  the  concentration  in  the  stagnant  volume  [A//L3],  C  is  the  concentration  in  the  mobile 
region,  A  is  the  fraction  representation  of  stagnant  to  mobile  [unitless],  and  K  is  the  constant 
describing  the  rate  of  exchange  between  stagnant  and  mobile  regions  [T-1].  Their  results  indicated 
the  total  mixing  observed  in  laboratory  experiments  was  the  result  of  more  than  just  the  dispersion 
mechanism.  (10:78) 

Van  Genuchten  and  VVierenga  continue  these  ideas  but  allow  for  a  sorbing  contaminant 
(36:473).  Using  their  notation,  the  following  expression  describes  contaminant  transport  within 
the  mobile  region  of  an  aquifer  in  cylindrical  coordinates: 

dCm(r,t)  _  A||V'(r)|  d~Cm(r,  l)  V(r)dCm(r,t)  9trnR,m  dC,ma(r,t) 
dt  Rm  dr 2  Rm  dr  0mRm  dt 

where  C;r„„(r,  t)  is  the  average  concentration  in  a  representative  volume  of  the  immobile  region 
[M / A3],  Ai  is  the  dispersivity  [L],  0m  and  #,m  are  the  porosities  of  the  mobile  and  immobile 
regions,  respectively  [unitless]. 

d'he  effect  of  th-se  immobile  zones  on  the  mobile  contaminant  concentrations  L  controlled 
by  the  ratios  of  the  porosities  of  the  two  zones  and  the  ratios  of  their  retardation  factors.  The 
retardation  factors  are  now 


Rm  —  1  + 


/  P  !<d 

0,n 


(3.7) 


for  the  mobile  zone  [unit less]  and 


Rim  —  1  + 


;i  -/)  p  k, 

0,m 


(2,8) 
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for  the  immobile  zone  [unitless].  The  term  /  denotes  the  fraction  of  sorption  sites  in  contact  with 
mobile  fluid  [unitless]. 

Rao  and  others  (29:684)  measured  diffusion  in  porous  ceramic  spheres  for  some  chemicals 
that  do  not  adsorb  in  an  attempt  to  describe  soh’te  transport  in  aggregated  porous  media  with 
distinct  mobile  and  immobile  zones.  A  Fickian  diffusion  mode)  was  used  to  fit  the  experimental 
data  The  results  pointed  favorably  toward  a  diffusion  model  to  describe  solute  transfer  between 
these  porous  spheres  and  the  mobile  region  (29:687). 

Nkedi-Kizza  and  others  examined  more  breakthrough  concentration  curves  from  laboratory 
experiments  to  explain  the  tailing  of  both  nonsorbing  and  sorbing  pollutants  (25:471).  They  coupled 
the  standard  advective/dispersive  model  presented  earlier  with  Fick’s  second  law  of  diffusion.  Their 
calculations  agreed  favorably  with  measurements  in  the  laboratory,  verifying  the  utility  of  assuming 
diffusive  solute  transfer  to  describe  transport  in  aggregated  porous  media  (25:475). 

Sudicky  and  others  (32)  put  this  model  to  the  test  in  a  laboratory  experiment  by  creating  an 
aquifer  of  thin  sand  sequestered  by  layers  of  silt.  Their  results  indicated  that  the  effect  of  molecular 
diffusion  within  the  silt  zones  was  more  important  than  the  dispersive  effects  in  the  mobile  zone  for 
understanding  the  fate  of  the  contaminant.  The  longitudinal  diffusion  (dispersion)  was  of  secondary 
importance. 

Miller  and  Weber  concent  rated  their  effort  on  the  rate-limited  effects  of  diffusion  by  examining 
a  variety  of  mathematical  models  for  predicting  contaminant  fate  and  transport  in  the  light  of  labo¬ 
ratory  experiments.  They  concluded  that  the  process  of  sorption  in  their  systems  was  rate-limited, 
and  the  best  fit  was  obtained  by  models  that  incorporated  intraparticle  diffusion  (24:243,2*50). 

The  equation  frequently  used  to  describe  Fickian  diffusion  within  the  immobile  zone  is 

dC;m(r,  z,t)  _  De  d 
dt  z‘'-lRimdz 

where  C\m  is  the  concentration  of  the  contaminant  in  the  immobile  zone. 

The  retardation  factor  R, ,n  accounts  for  contaminant  sorption  within  the  immobile  zone,  and 
is  called  the  immobile  zone  retardation  factor.  The  diffusion  coefficient  in  tlw  i,,miot>iie  zon,  i.-  I), 
[A'//’].  The  term  u  defines  the  type  of  immobile  region  geometry: 


,v-\ 


dC',m{r.  :,<)! 


8z 


(2.9) 
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i/  =  1  for  Layered  Immobile  Region  Geometry 
r  —  2  for  Cylindrical  Immobile  Region  Geometry 
v  —  3  for  Spherical  Immobile  Region  Geometry 

To  couple  this  equation  with  the  mobile  region  mass  balance  equation  (Equation  2.0).  volume' 
averaging  is  used.  The  amount  of  contaminant  entering/leaving  the  mobile  region  must  equal  the 
amount  of  contaminant,  leaving/entering  the  immobile  region.  The  amount  of  contaminant  in  tin? 
immobile  zone  is  given  by 

C,ma(rJ)=£  f  :'-'flm(r,.-,t)rf:  (2.10) 

&  Jo 

where  b  is  the  halfwidth  of  a  typical  immobile  zone  [/,],  and  u  defines  the  type  of  immobile  region 
geometry.  One  of  the  findings  of  Nkedi-Kizza  and  others  (25:475)  was  that  a  reasonable  range  of 
sizes  and  shapes  of  aggregates  could  be  approximated  by  a  representative  spherical  aggregate. 


2.4  The  .\fodel  Equation  Set 

The  equation  for  the  mass  balance  (Equation  2.6)  needs  the  knowledge  of  what  transpires  in 
the  immobile  zones,  to  determine  the  last  term  on  the  right  hand  side  of  the  equation. 

The  governing  equation  for  the  immobile  zone  depends  upon  whether  u  is  one,  two,  or  three 
For  v  —  1  Equation  2.9  becomes 

dCim(r,zJ)  De  d2Ci,n(r,  :,t) 
m  /?„„  d~J- 

and  tin'  volume  averaging  equation  (Equation  2.10)  becomes 
1  fb 

C'ima(r.t)  =  -  J  C'„n(r.z.t)dz  (2.12,1 

These  equations  (2.(3,  2.11,  and  2.12)  form  a  complete  set.  describing  advection  and  dispersion 
in  a  radially  symmetric  aquifer  with  immobile  layers  of  clay  or  other  layered  sorbate  (see  Figure  2. 1 ). 


When  v  —  li  Equation  2.9  becomes 


Figure  2.1.  Layered  Immobile  Zones  in  an  Aquifer 


and  Equation  2.10  becomes 
3  rb 

(  l ma(n.  11  ~~  J  ~  \  0  dr 


( 


These  equations  (2.0,  2.13,  and  2.M)  describe  spherical  immobile  regions  in  the  aquifer  u 
the  contaminant  may  undergo  adsorption  (see  Figure  2.2). 

These  equations  form  the  basis  of  the  numerical  model  examined  in  the  next  chapter. 
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Figure  2.2.  Spherical  Immobile  Zones  in  an  Aquifer 
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III.  Numerical  Methods 


3.1  Problem  Statement 

la  the  last  chapter  we  looked  at  the  theoretical  and  experimental  evidence  for  including;  rate- 
limited  sorption/desorption  physics  in  our  contaminant  transport  model  This  chapter  pro*  nts 
the  numerical  implementation  of  these  equations.  This  implementation  uses  numerical  quadrature 
techniques,  the  Finite  Element  Method  (FEM)  and  Finite  Difference  Method  (FDM)  to  transform 
the  partial  differential  equations  into  an  algorithm  suitable  for  t hi’  Fortran  language. 

■  I.J  Program  Design 

Beginning  with  Equations  2.9.  2.10.  and  2.(5  as  representing  either  the  layered  or  sphere  d 
geometry,  the  preferred  solution  would  be  a  transform  that  combines  the  two  regions  and  tin- 
three  equal  ions  into  a  single  equation.  The  Laplace  t  ransform  was  attempted,  but  the  rest  re  m  a- 
of  a  formula  representation  of  the  levels  of  contaminant  in  both  regions  was  uuacceptahh  \V,’L 
this  restriction,  the  model  could  not  simulate  pulse  pumping  These  equations  are  solved  for  tie 
special  case  of  a  constant  level  of  the  contaminant  from  the  well  center  out  an  arbitrary  distance 
and  then  down  to  zero  (lit).  This  Laplace  transform  solution  will  be  discussed  in  Chapter  1 

The  simplest  and  most  obvious  solution  was  to  keep  arrays  of  gridded  data  for  each  ivge  n 
Figure  3. 1  shows  the  representation  implemented  in  this  thesis.  The  mathematical  formula  h  r 
tin'  seepage  velocity  in  the  mobile  region  presents  a  fictitious  singularity  at  the  well  cent .  r  E< 

overcome  this  problem,  the  grid  begins  a  slight  distance  from  the  well  center:  an  approxima’  e  n 

to  I  lie  boundary  of  the  well  itself.  The  grid  for  the  immobile  zone  is  more  elaborate  than  lh.it 
for  the  mobile  zone  Inv.msc  of  the  need  to  maintain  a  sufficient  number  of  grid  points  to  pr.  q  -  rU 
represent  a  gradient  profile  for  the  Fickian  diffusion. 

The  model  will  step  forward  within  the  immobile  region  a  small  amount  of  time  ( s. .  1  v ; :  i  e 
Equation  2.9.  and  then  integrate  (solving  Equation  2.10).  These  new  values  will  then  be  u-  I  t 

calculate  new  values  iu  the  mobile  zone  by  stepping  forward  in  time  using  Equation  2  e,  1  in- 

pro-ess  will  be  repea'-d  eetil  the  final  time  ts  reached. 
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i.  J.  I  Data  Flow  Oriented  Design  In  order  to  implement  pulse  pumping.  tin-  model  i ••  1.- 
t lie  ability  to  stop  the  flow  rate  after  a  specified  t  hue  ami  begin  again  with  reduced  or  increas. •> i  If  av 
rates,  and  perhaps  adjust  other  parameters  like  the  duration  of  the  simulated  execution.  Halle  r 
than  have  the  model  run  interactively  and  solicit  information  periodically  from  the  user,  it  will 
create  an  output  file  that  contained  the  same  information  as  the  input  file.  With  this  design,  tie 
model  could  cycle  on  its  own  data,  keeping  the  output  contaminant  levels  for  both  regions,  grid 
si/es.  and  other  grid  information  from  one  execution  and  using  this  as  input  to  the  next  run 

Because  of  the  large  number  of  input  parameters,  two  input  files  were  chosen:  one  for  the 
parameters  that  were  subject  to  change  and  tuning',  and  one  for  the  grid  and  contaminant  |e\cl 
informal  ion 

I  hese  files  are  shown  in  the  data  flow  diagram  (Figure  3.2)  and  the  formats  of  these  files  are 
slu  iwn  in  the  document  at  i<  >n  leader  for  the  r<  >ut  dies  that  read  them  (( i  HARMS  and  ( !  FT  IK  'Mi 

t  ~  J  I  ransfoi  ni/'l  nm  on  turn  Analysis  transform  analysis  was  chosen  as  an  overall  design 
methodology  because  of  the  nature  of  the  problem  .liferent  (low  from  the  input  files,  m  iv  lube- 
computed,  and  c  fie  rent  flow  for  the  data  to  the  output  file  ('2*v2Su ) 

I  lie  main  routine  receives  all  tin-  input  information  from  (ISINK).  passes  all  this  to  I'lll  HI 
wlmh  returns  the  final  data  (see  Figure  3  3  for  the  overall  implementation  based  on  transform 
analysis)  I  lie  mam  rout  me  then  passes  the  final  data  out  to  SAYl  IK '  which  creates  the  output 

file 

Ivys  m  the  input  data  deter  nine  the  geone-i  ry  of  the  immobile  zone,  and  (lie  choice  of  fin  n . 
differences  or  finite  elements  in  l lie  mobile  zo ne  I  hese  possibilities  (jiiite  naturally  fall  out  as 
transactions  (2K:27.">).  Figure  ,'l  t  shows  the  transact  ion  analysis,  with  I’RKFl  as  the  transaction 
center  and  several  choices  beneath  it  for  the  differ,  nt  possibilities  mentioned  above. 


MODKL1T 


Figure  3.3.  MODKL1T  Transform  Analysis 


Figure  3.1  MODFUT  Transaction  Analysis 


:i.:i  Numerical  Methods 


The  following  sections  describe  the  procedures  employed  for  mapping  the  partial  different  ia! 
and  integral  equations  into  numerical  algorithms. 

3.3.1  7’Af  Immobile  Zone  The  immobile  zone  grid  (shown  in  Figure  3.1)  is  somewhat  de¬ 
ceptive:  the  grid  spacing  in  the  radial  direction  is  on  the  order  of  meters,  while  the  grid  spacing 
ill  the  z  direction  is  on  the  order  of  millimeters.  In  other  words,  each  column  in  the  immobih 
zone  grid  is  virtually  isolated  from  every  other  column.  The  contaminant  diffuses  in  a  column,  but 
not  between  columns.  This  should  be  intuitive  because  each  “column”  is  actually  a  representative 
layered  or  spherical  aggregate  of  the  immobile  material. 

The  immobile  zone  could  consist  of  spherical  aggregates  or  layers.  A  column  in  the  immobile 
zone  represents  the  concentrations  in  a  halfwidth  of  a  layer  or  the  concentration  profile  along  the 
radius  of  a  typical  spherical  aggregate,  depending  on  whether  we  are  studying  layered  or  spherical 
zones,  respectively. 

3.3  1.1  Diffusion  in  Layers  The  derivative  on  the  right  hand  side  of  Equation  2.11  is 
replaced  by  a  numerical  approximation  based  on  the  Taylor’s  series  expansion  about  the  point  it 
for  an  arbitrary  analytic  function  /(te): 


/(<<■  +  Ate)  =  f(w)  + 
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and 
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These  are  added  and  the  resulting  infinite  series  truncated  to  produce  the  convergent  uum.  nca! 
.approximation  for  the  second  derivative  with  gridded  data 


m,  ,  /(te- Ate)  +  2/(  te)  + /( te  +  Ate)  , 

/  (te)  =  - - - r— 5 -  +  t>( Ate  ) 

At  e- 

Wnh  this  approximation  applied  to  - - "r‘  liquation  2.11  becomes 
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dCim(r,z,t) 

dt 


Dc  C,m(r,  :  -  Ar,  t)  +  2C„n(r,  s,t)  +  Clm(r,  ;  +  A zj) 

R„n 


Eliminating  the  derivative  in  time  by  following  the  method  of  Crank  and  Nicolson  (.'$8:107'))  by 
integrating  from  t  to  t  +  A t: 


f+*‘  dCim(r,  :,t) 

Jt 


De  C\m(r,  z  -  A r)  -1-  2C,m(r,  z,  t)  +  Clm(r,  -  +  A;,  r) 
IU,n  A:3 


dr 


and  using  Tie  trapezoidal  quadrature  ride  to  get  C,m  at  the  future  time: 


C„u(r.z,l  +  A/)  =  C,m(r,;,t) 


[C  im(r,  :  —  Ac,  t  +  At )  +  2C,m(i\  /  +  At )  +  C'i,,,  (r,  z  +  Ac ,  I  +  A  / ) 

+  Gm(t  -  —  A:,  f )  +  2Cjm(r,  1. 1)  +  Cim(r,  :  +  A;,  1)]  (■$.(>) 


'l'liis  method  also  requires  boundary  conditions  for  the  surfaces  of  the  immobile  layers  at  t 
and  t  +  At.  The  values  at  the  boundaries  of  these  layers  are  the  contaminant  concentrations  in  the 
mobile  region. 

I  bis  implicit  method  has  the  following  benefits  (17:132.152):  it  is  absolutely  stable,  it  does 
not  change  amplitude  of  the  physical  mode,  and  has  no  computational  mode!  Its  drawbacks  are  i  lie 
need  to  do  a  matrix  inversion  and  it  retards  the  phase  slightly.  The  routines  IMA'EIH,  IM.VE3H. 
IMA  I  If,  and  IMAT1M  all  create  matrices  using  this  method,  and  invert  these  matrices  once:  no 
further  matrix  inversions  are  needed  during  the  execution. 

I  he  grid  points  are  numbered  with  point  1  nearest  the  mobile  zone,  and  increasing  towards 
I  lie  center  of  the  immobile  zone.  1  he  contaminant  concentrations  at  one  side  of  the  imtlcbih 
zone  match  those  at  the  other.  Ellis  is  because  the  concentration  in  the  mobile  zone  has  no 


vertical  component.  For  this  reason,  the  concentration  profile  of  the  contaminant  in  the  immobile 
layers  has  a  symmetry  around  the  renter.  Equation  3.(5  need  only  be  expanded  for  the  halfwidth. 


with  a  reflection  across  the  center  (at  z  =  0).  Taking  z  =  (i  -  l)Ac  for  i  =  1,2,  ,  A/  and 

Cim,(r.  0  =  Cim(r,  (»  —  IjAz,  f),  the  matrix  form  of  Equation  3  6  is: 

B\ Cm(r,  t  +  At)  =  BiCtm(r,t)  +  U^r.t)  (3  7) 

the  matrices  and  vectors  are  defined  as: 


B 

fli.-,.. 


1  + 


D'£t 

R,m£r2 


-D'&t 

2R,m£r2 


-D'&t 
2/?,m  Ar2 


•  i,=  1,2,  •  ,M-  1 

«,=  2,3,  .,Af 
i,=  2,3,  -,.Vf 


(3  S) 
(3.9) 
(3.10) 


fla... 

So..,., 


_  De£t 
Ar2 

Dt£t 
2/^m  Ar2 
D;  At 
2/2,m  Ar2 


i,  —  1,2,  -,M  -  1 

»,=  2,3,  ,  A/ 

»,  =  2, 3,  •  ,  Af 


1 


De  At 
2. Rim  Ar2 


[Cm(r,t)  +  Cm(r,t  +  At)] 


i  =  1 


(3-11) 
(3  12) 
(3  13) 

(3- Id) 


[  0  otherwise 

The  boundary  conditions  mentioned  above  are  incorporated  into  the  V'i  vector.  This  scheme  requires 
the  knowledge  of  Cm(r,  t  •+  At)  which  is  not  available  until  after  several  more  procedural  steps.  An 
approximation  to  Cm(r,  t  +  At)  was  made  using  the  available  data  and  linear  extrapolation: 


Cm(r,  f  +  At)  =  2Cm(r,  f)  -  Cm(r,  f  -  At) 

These  equations  are  implemented  in  routine  IMAT1B  in  the  program. 


(3  15) 


3.3.1.S  Diffusion  «n  Spheres  The  form  of  the  Fickian  diffusion  equation  for  spheres 
differs  slightly  from  that  for  layers.  The  additional  term  in  Equation  2.13  is: 

-  dc>m(r.  z.t)  (3  16) 

z  dz 

Using  the  Taylor’s  series  expansions  above  (Equations  3.1  and  3.2),  and  subtracting  instead 
of  adding,  to  obtain  the  approximation  for  the  first  derivative: 
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=  /(u_4- J|u.  ^u-)  +  ^ 

2Au> 


This  is  included  in  Equation  3.4  resulting  in: 


r'+A(  dCtm(r,z,T) 


ft+*‘  De_  Cim(r,  z  -  Az,  r)  +  2Cim(r,  c,  r)  +  Cim(r,  c  +  A--,  r) 

J,  Run  Ac7 

1  Cim(r,  c  +  A:,T)-Clm(r,:-  A:,t)]  j 

+  - - - - - - T- - - 

2  A  2 


which  is  integrated  using  the  trapezoidal  quadrature  rule  to  obtain  the  numerical  scheme 


Cim(r,:,t  +  At)  =  Cjm(r,  2, 0  +  Af 


2ftraA:2 


[C<m(r,  2  -  Ac,  <  +  At)  +  2Cim(r,  2,  t  +  A<)  +  Cjm(r,  2  +  Ac,  /  +  AO 
+  -  (CVm(r,  2  +  Ac,  t  +  At)-  Clm(r,  2  -  Ac,  t  +  At)) 


+  Om(r.  2  —  Ac,  0  +  2C’im(r,  c,0  +  Ci,n(r ,  c  +  Ac,  /) 


+  -  (£'»»('•.  2  +  Ac,  0  -  Cim(r.  2  -  A2.0) 


In  vector  form  this  is  identical  to  Equation  3.7  with  slightly  different  constituents  for  the 


matrices. 


Si...-,  = 


V  Rimh'J' 

-DeAt  fi  -  2 
2  /?,•„,/, 2  \  *  —  1 
-D,A<  /  j 
2R„nh-  1,-1 


j  =  2,3,  -  -  -,  A/ 


=  2, 3,  •  •  • .  A/  -  1 


i  =  2, 3,  •  •  • ,  A/  -  1 


(3.20) 


(3.21) 


(3.22) 


V  KimAV 

DeA<  -  2 

2 2  V^T 

DfAt  /  i 
2Rimh-  \i  -  1 


/  =  2,  3,  ,  ,U  -  1 


,  =  2,3,  -,  M  -  1 


=  2,3,  ■  ■  • ,  A/  -  1 


(3  2  1) 
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Vi.(r,/)=  { 


(3.2b) 


EE  '  [Cm(r,t)  +  Cm{r,t  +  A/.)]  <=1 

^  m  ^  ~  “ 

^  0  otherwise 

When  i  —  1  these  formulas  can  not  he  computed,  but  the  first  grid  point  in  the  mobile  region 
begins  .  l^r  from  the  well,  yielding  a  slightly  different  value  for  the  B\2  ,,  B j,  B-2l  ,,  and  B-,i  , 
components.  This  is  implemented  in  routine  IMAT3B  in  the  program. 

3.3.2  The  Volumetric  Averaging  To  get  the  time-rate-of-change  of  the  pollutant  concen¬ 
tration  in  the  immobile  zone  we  integrate  over  the  halfwidth  of  the  layer  or  over  the  radius  of 
the  sphere.  This  change  is  what  makes  its  way  into  the  mass  balance  equation  {liquation  2.1)) 
to  account  for  the  chemicals  entering  and  leaving  the  mobile  zone.  The  equation  governing  tin.’ 
volumetric  averaging  (Equation  2.10)  is  approximated  for  each  of  the  two  types  of  immobile  zones 
by  the  quadrature  formula  known  as  Simpson's  rule  (3S:369). 


3.3.2. 1  Volumetric  Averaging  in  Lagers  Equation  2.12  gives  the  volumetric  average  of 
the  pollutant  in  a  typical  layer  of  the  immobile  zone.  The  layered  medium  provides  t  he  simplest  form 
of  Equation  2.10  and  the  conversion  of  this  integral  to  a  numerical  quadrature  is  straight  forward. 
Simpson’s  rule  requires  an  odd  number  of  points.  For  example,  given  a  funct  ion  /  which  is  integralde 
on  the  interval  [a,  6]  : 


I 

I 

l 


f  fW 

J  a 

(b  -  a) 

3 

(h  -  a) 


elw 


f(ei)  +  4/  (  — - —  ]  +  f(b) 


(5 


/(a)  +  4/  a  + 


(b  -  a) 
4 


+  2/  a  + 


2(6  -  a) 


4 


+ 


4/  (  a  +  —  4  —  )  +  f(b) 


(3.27) 

(3.2S) 


(3.20) 


But  with  symmetry  about  the  origin,  the  number  of  grid  points  in  a  column  in  the  immobile  zone 
can  be  even  or  odd.  The  point,  at  ;  =  0  is  the  renter  point,  in  the  region,  with  grid  points  extending 
in  both  directions,  necessitating  an  odd  number  of  points  total.  The  algorithm  implement'  d  in 
.\  It 'IMA  makes  use  of  this  argument  to  integrate  over  the  halfwidth. 
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3. 3. 2. 2  Volumetric  Averaging  in  Spheres  The  argument  for  spheres  is  almost  the  same 
as  that  for  layers,  except  that  here  the  number  of  grid  points  requested  in  the  immobile  region  must 
be  odd.  The  model  checks  to  make  sure  the  requested  number  of  grid  points  will  work  with  this 
algorithm  in  the  GSINFO  routine. 

The  integrand  of  Equation  2.14  contains  a  z 2  term  which  complicates  matters  slightly.  The 
center-most  grid  point  of  a  sphere  doesn’t  contribute  to  the  total  volumetric  average,  as  is  shown 
in  this  example  formula  where  A:  =  |  (i.e.,  5  grid  points  in  the  immobile  zone): 

1  =  fr  [/(*)  +  4/( 3A-~)  +  2/(2Az)  +  4/( Az)  +  /( 0)]  (3.30) 

where  f(z)  =  z2Ci,„(r,  z,  t)  for  fixed  r  and  /,  and  the  result  is 

A  , 

I  =  ~  [fc2Cim(r,&,0  +  4(3Az)2Clm(r,3Az,f)  +  2(2Az)2Cim(r,2Az,l)+ 

4(A.-)2C,m(r,Ar,<)+(OAz)2C,m(r,0,<)]  (3.31) 

3.3.3  The  Mobile  Zone  Mass  Balance  Equation  The  mass  balanced  equation  (Equation  2.0) 
is  solved  numerically  by  two  methods,  the  finite  element  method  and  finite  differences.  The  finite 
differences  are  easier  to  understand,  and  were  implemented  as  a  check  against  the  finite  elements. 
Comparisons  of  these  two  methods  with  the  same  input  data  are  made  in  Chapter  4. 

3.3.3. 1  Finite  Differences  in  the  Mobile  Region  The  procedures  for  the  finite  difference 
method  involve  replacing  the  deri vales  in  the  equation  with  approximations  to  the  derivatives  based 
on  Taylor’s  series  expansions.  The  partial  derivatives  with  respect  to  r  of  Equation  2  0  are  repln'S-il 


d2Cm(r,t) 

dr2 

dC'mjrJ) 

dr 


C m  ( r  —  A r,  t )  +  2 C  m(r,  t )  +  Cm(r  +  A7\  t ) 


Ar2 

C'm(r  +  Ar,t)  -  Cm(r  -  A r.t) 
2A  r 


resulting  in 


dCm(r,t)  dim  Rim  dCjmafr,  t ) 

dt  ‘  0mRm  dt 


(3.33) 
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At\V(r)\  I- Cm(r  -  Ar,  t)  +  2 Cm{r,  t)  +  Cm[r  +  Ar,  t) 
Rm  [  Ar2 

_  V(r)  r C'm(r  +  A r,t)  -  C'm(r  -  AM)'  ^  . 
Am  2Ar 


(3-3-1) 


Integrating  this  from  t  to  t  +  At  and  rearranging  yields: 


C'm(r,t+At)  =  Cm(r,t) 


'fa  ir  [Cimairj  +  At)  -  C,ma(r,t)] 

LAt  Um  ilm 

A^  [^/|V(r)l 
2  ^ 

f  C m(r  —  &r,t  +  At)  +  2Cm(r ,  t  +  At )  +  Cm (r  At-,  t  -f  At ) 

\  Ar2 

V(r)  f  Cm(r  +  Ar,  t  +  At)  —  Cm(r  —  Ar,  t  +  At )  \ 

Rm  V  2Ar  ) 

A,|V'(r)l  /C,„(r-  Ar,Q  +  2Cm(r.<)  +  Cm(r  +  Ar,t)\ 

V  Ar2  ) 

V(r)  /  Cm(r  +  A  r,t)  -  Cm(r  -  Ar,  t)Y 
Am  V  2Ar  J 


Taking  r  =  (i  -  l)Ar  for  i  —  1,2,  ,  At  lr>  tt>;«  equation  gives  the  vector  equation: 

Pi  Cdt  +  At)=  P2Cm(t)  -  j-t  [C,mo(t  +  At)  -  C,ma(t)1  (3.3(5) 

where  /i  ■-  ~ and  Cm,(t)  =  Cm([»  -  l]Ar.f)  for  i  =  1,2,--,  At. 
f'm 

The  values  for  Clma(r,  t)  are  available  at  t  +  At  because  they  depend  upon  C\m(r,  :,l)  only, 
and  the  algorithm  makes  the  linear  extrapolation  to  approximate  Cm(r,  t  +  At).  The  tri-diagonal 
matrices  P\  and  P->  are  initialized  in  routine  IMAT1F,  and  individual  components  of  these  matrices 
are  defined  here: 


/  =  2,3,... .At-  1 

(t  -  1) 

■— —  +  —  i  =  2,  3,  ■  •  ■ ,  A'  —  1 

(i-1)  At 


(-C+Q 

(i-D 


i  =  2,  3,  ,  At  -  1 


(3.3!)) 


(C  +  Q 

(«-  1) 


i  =  2.3,  ,  At  -  1 


(3  10) 


And  for  i  6  {■'!.  1 .  ■  ,  A" } : 


Mr) 


-  —  (j  —  2)  (»-  2)h  <  r  <  (i-  l)h 

h 

'  -7-  +  *  (i  -  l)/i  <  r  <  (i')A 

h 

0  otherwise 


0>+i  (,-)  —  ^  — — I-  (i  +  1 ) 

h 


( i  —  \  )h  <  r  <  ( i)h 
(i)h  <  r  <  (j  +  \)h 
otherwise 


<Pi-i(r)  =  { 


[-U-3) 


o 


(i  -  3)h  <  r  <  ( i  -  2 )h 
( i  -  2  )h  <r<  (i  -  1  )h 
otherwise 


with  derivatives 


o\  (r)  =  < 


0.9/i 

0 


(0.1  )/i  <  r  <  h 
otherwise 


o',(r)  =  { 


1 

0M 

-A 

h 


0 


(0.1  )/i  <  r  <  h 
h  <  r  <  2h 
otherwise 


O'ii  r) 


(i  —  2)/i  <  r  <  (i  -  l)/i 
(»  —  \  )h  <  r  <  (i)h 
otherwise. 


'■  +  1 


>•) 


(i  -  \)h  <  r  <  (i)/i 
(i)/i  <  r  <  (j  +  l)/l 


0  otherwise 


1 

h 


(i  -  3)/i  <  r  <  (i  -  2)/i 
(  /  -  2)/i  <  r  <  (J  -  1  )/i 


(3,lo) 


(O-IO) 


(3,17) 


(3,1S) 


(3.19) 


(3.30) 


(3.31) 


(3,32) 
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0  otherwise 


Taking  Equation  2.6  and  multiplying  by  <t>j(r)  for  each  j  =  1,2,  -  A'  and  integrating  with 
respect  to  r  over  the  domain  [0 , R]  yields: 


t 


dC,n(r,  t) 


a,  «r)Jr  =  7 it  l  |v’wl 

-  if 
L 


dr 2 


4>j{r)dr 


,dCm(r,t) 

V  (r) - — - 0j(r)dr 

—  3  I  —  im3^  — - <t>j(r)dv  for  i  ~  1,2,  •••,A' 


dt 


Integration  by  parts  of  the  second  derivative  term  produces: 


(da.i) 


/V(’’)l 

*/0 


a2^(r,o ^  /  w  wacm(r,oJ  ,  ,  * 

0j(r)dr  =  (V(r)| - - - <2>;(r) 

or  o 


9r- 


-  L 

-l 


RdCm(r,t) 

- <£j(r)|V  (r)|dr 


0Cm(r,t) 

Or 


<f>j(r)\V'(r)\(tr 


(d. ad) 


Applying  the  definition  of  the  <?j(r)  eliminates  the  evaluation  term.  Substituting  this  result  back 
into  Equation  3.5d  yields: 

A,  dCm(r,t ) 


rH  dCm(r,t)  .  ,  At  ffidcm(r,l)„,  . . 

L  =  -s:/0 


- 

R,n  Jo 

iC  l 

-  >fl- 

Jo 


dr 

KdCm(r.t) 
dr 


<Mr)|V"(r)|rfr 


i  •/  i  oc^ ( t*.  I)  w 
l  (;•) - — - Oj(r)dr 


Or 


OC 

ima  (r.  0 
dt 


<t>j(r)dr  for  j  =  1 , 2,  ■  ■ ,  A’ 


(d.  •’>.’>) 


The  variables  Cm(r,l)  in  this  equation  are  replaced  by  the  chapeau  functions  shown  above 
N 


Cm{r,t)  =  ^  r‘« (f  )<?«('•) 


(d.ati) 


i  =  1 


giving: 


v 

/  (r)Oj(r)dr  = 

./()  , 


J  /■  ft  ^  ^ 

-77^-  /  V«,(Mo:(r)0;(r)|r(r)|dr 

“m  Jo  |=] 

*  m  Jo  _  , 


rli  -V 

./o  ,-l 


)0;(r)V'(r*)f/r 


-i 


« 


,  I  SGmu(r.<)  ,  ,, 

i  /  - — - Oj(r)dr 


(•'i  m! 


for  j  —  1 . 2,  •  •  • ,  A’.  This  produces: 


A 


/  +  3  [ 

Jq 


dt 


0(r)dr  =  -  —-A2  0-  At  a  -  A,  a 

‘^773  *»m  *  i  rn 


(d  AS  J 


where  o(t)  =  [o ,  ( / ) ,  cv2(<) . or,v(/)]'r  and  «'(0  =  [«'(/),  a'2(<) _ a'A,(/)]7\  and  the  basis  func- 

t'oti  vector  o(r)  =  [d>i  (r ),  <p?(r),  -  •  • ,  <3,\(r)]r ,  with  the  matrix  components  for  the  four  matrices 
A .  A_>,  At.  and  At  defined: 


A„, 

=  f  0,(r)dj(r)dr  i  =  1, 2,  •  • 

J  () 

- ,  A’  j-  1,2,  -.A' 

(:i .  A : ) ) 

r'< 

A-. , 

=  /  o;(r)o;.(r)|V(r)|rfr  »  =  1.2. 

Ju 

.  A'  j  =  1 . 2 ,  ■  ,  A' 

(d  I'.m 

fl< 

A„, 

=  /  <P,,('-!oJ('')|T'(e)|r/r  /  =  1.2. 

Jo 

■  ■ ,  A’  j  -  1 , 2,  .  A' 

(d.i'd.i 

-  /  e>!(r)e>j(r)V'(r)r/r  i  =  1,2. 

J{) 

.  .V  j  =  1 , 2 ,  ■  ■  • ,  .V 

(d  02) 

rite 

method  of  Crank-Nicolson  (trapezoidal  cpiadrature  for  the  time  derivative) 

is  Used  l.  , 

remove  the  integration  in  time  to  obtain: 

<*(/  +  At) 


M.  A  A.i 

b  rrrr~  (Aj  +  At)  +  — 


At  2IL 


2/f„ 


A  _  -dr 

At  2/?„, 


(  A_>  T  At )  —  - 


A 

2ll„ 


£l(0 


-  ,i 


["  (r.t) 

Jo  dt 


o(  r)dr 


which  is  further  simplified  by  defining 


O, 


(dr.  1 1 


o. 


A  A  A( 


i  :t. 


yielding: 


d-  1 5 


Q\<±(t  +  Af)  =  Q2o(0  -  jj  ■-C'm.^--6{r)dr 

Those  matrices  are  initialized  in  routine  IMAT1M  in  t lie  program  MODKLIT.  (s< 
this  Thesis). 


i  :i 

■e  Part 


1  (i 


IV.  Results 


4-1  Introduction 

This  chapter  presents  an  evaluation  of  the  computer  model  developed  in  Chapter  it  Su.. , 
there  exists  a  paucity  of  suitable  field  data  for  comparison,  comparison  will  be  made  with  !:■ 
transform  solution  provided  by  Dr,  Goltz.  An  internal  comparison  will  also  be  present.  |  t,, 
demonstrate  the  relative  effect  of  the  advection,  dispersion,  and  diffusion  mechanisms. 

The  results  of  several  executions  are  presented  here.  The  model  execution  produces  eteni.  i~ 
quantities  of  data  that  almost  require  graphical  depictions  to  make  them  understan  lalT  I  !■■ 
following  sections  include  many  such  pictures  showing  the  progress  and  fate  of  the  contaminant  In 
many  instances,  the  output  from  one  execution  is  fed  back  into  the  next  execution  to  demon-trite 
the  feasibility  of  simulating  pulse  pumping. 

.{■J  Comparison  With  the  Analytic  Solution 

The  computer  model  solves  the  same  equation  set  list'd  by  <!o!tz.  but  surpasses  the  abtli'ti*  - 
of  his  solution  by  allowing  arbitrary  initial  conditions,  and  pulse  pumping.  The  cost  of  thi-  ad  !•  1 
flexibility  is  the  accuracy  of  the  solution  as  time  progresses.  With  ident  ical  initial  tlata  both  m.  ■  !■  I- 
should  produce  similar  results.  The  analyt  ir  solution  solves  t  he  problem  for  all  t  imes  si  mu  It  an..  ■  > :  - 1  \ 
and  tlie  numerical  solution  steps  forward  in  time,  and  is.  therefore,  subject  to  compoim  ling  <  n  r~ 

I  Laytnil  Immobile  Mahrinl  The  results  of  one  comparison  are  shown  m  I .  i  i  ■ !  •  !  1 
wlnre  the  column  labeled  TKM"  is  the  finite  element  solution,  the  coin  tut:  labeled  T  MM  i>  (me 
dillerences.  and  the  "Anal"  is  the  analytic  solution  provided  by  Goltz  (Kit.  1  hose  vain-  -  ar.  tie 
concentration  of  the  pollutant  at  the  well. 

I  lie  initial  cont  a  min  ant  concent  rat  ion  distribution  was  a  const  ant  1 .0  out  2*  meter.-  If.  m  t  !;■ 
well  center  and  then  downward  to  a  value  of  0.0  for  hot h  mobile  and  immobile  zones  1  he  mi m.  :  I 
zt ate  type  was  ‘layered’,  the  other  input  data  is  shown  in  (able  1.2. 

The  TIM  solution  is  reasonably  close  to  the  analytic  solution,  and  in  fact.  onl\  de\r,;.~ 
by  n  TV  A  at  100  days!  The  “l  KM  solution  more  closely  n'sembled  the  analytic  solution  than  tie 


bl 
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Figure  4.1.  Mobile  Zone  Concentration  Progression,  Days  0-10,  Layered  Immobile  Material 

"FDM”  solution.  Part  of  the  deviations  between  these  solutions  is  caused  by  the  initial  profile  of 
the  concentration.  Whereas  the  analytic  solution  drops  the  concentration  from  1.0  to  0.0  instantly 
at  r  =  R,  the  numerical  solution  drops  it  linearly  between  two  grid  points.  This  makes  it  difficult 
to  match  the  initial  conditions,  especially  when  the  grid  spacing  is  over  2  meters,  as  in  the  first 
test. 

Figures  4.1  through  4.6  show  the  level  of  the  contaminant  in  the  mobile  region  as  time  pro¬ 
gresses.  Each  of  these  figures  represents  10  days  time,  with  the  beginning  profile  in  the  foreground 
and  ending  in  the  background.  The  magnitude  of  the  contaminant  is  shown  by  the  axis  of  the 
ordinate  and  the  radial  distance  out  from  the  well  is  shown  along  the  axis  of  the  abscissa.  The 
greatest  and  least  levels  of  the  contaminant  shown  on  the  graph  are  printed  on  the  right  rear  of 
the  projection.  The  angle  of  the  projection  is  the  angle  the  perspective  appears  to  make  relative 
to  the  abscissa.  The  numbers  along  the  abscissa  are  the  number  of  grid  points  in  the  radial  dimen¬ 
sion,  with  1  being  the  point  at  the  well  boundary.  The  numbers  along  the  right  side  of  the  graph, 
following  the  perspective,  are  the  number  of  iterations  displayed  (not  the  number  of  iterations  for 
the  calculations). 

As  these  figures  demonstrate,  the  contaminant  is  drawn  toward  the  well  and  extracted,  until 
the  level  drops  below  about  10%.  From  this  point  onward,  the  level  drops  much  more  slowly, 
until  it  appears  to  level  off  at  around  9.070.  Actually,  the  level  is  dropping,  but  the  contaminant 
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Figure  4.2.  Mobile  Zone  Concentration  Progression,  Days  10-20,  Layered  Immobile  Material 


Figure  4.3.  Mobile  Zone  Concentration  Progression,  Days  20-30,  Layered  Immobile  Material 


2 . 7E-01 
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Figure  4.5.  Mobile  Zone  Concentration  Progression,  Days  40-50,  Layered  Immobile  Material 


Figure  4.6.  Mobile  Zone  Concentration  Progression,  Days  90-100,  Layered  Immobile  Material 
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Figure  4.7.  Immobile  Zone  Concentration  Profile,  Day  10,  Layered  Immobile  Material 


Figure  4.8.  Immobile  Zone  Concentration  Profile,  Day  20,  Layered  Immobile  Material 

that  desorbs  from  the  immobile  regions  feeds  the  mobile  zone.  This  can  be  seen  more  readily  by 
examining  Figures  4.7  through  4.12  which  display  the  concentration  profiles  in  the  immobile  zone 
at  the  end  of  each  period  indicated. 

The  level  of  pollutant  in  the  mobile  zone  is  displayed  in  the  foreground,  with  the  profile  for 
the  immobile  zone  behind  it.  As  the  level  of  the  pollutant  in  the  mobile  region  decreases,  the 
gradient  between  these  two  regions  increases,  as  does  the  flux  of  chemicals  from  the  immobile  zone. 
At  the  end  of  this  100  day  test  (see  Figure  4.12),  the  levels  in  the  immobile  regions  are  still  quite 
high,  and  indicate  that  only  a  small  portion  of  the  total  contaminant  was  removed. 


Figure  4.9.  Immobile  Zone  Concentration  Profile,  Day  30,  Layered  Immobile  Material 


Figure  4.10.  Immobile  Zone  Concentration  Profile,  Day  40,  Layered  Immobile  Material 


Figure  4.11.  Immobile  Zone  Concentration  Profile,  Day  50,  Layered  Immobile  Material 


i 
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Figure  4.12.  Immobile  Zone  Concentration  Profile,  Day  100,  Layered  Immobile  Material 


Table  4.3.  Spherical  Immobile  Material  Comparison  Test 


Day 

Anal 

FEM 

FDM 

Day 

Anal 

FEM 

FDM 

0 

1.000 

1.000 

1.000 

30 

.0381 

.0518 

.0546 

2 

1.000 

1.000 

.9999 

40 

.0250 

.0359 

.0377 

4 

1.000 

1.000 

.9813 

50 

.0170 

.0258 

.0271 

6 

.8174 

.8187 

.8343 

60 

.0116 

.0188 

.0197 

8 

.4468 

.4493 

.5405 

70 

.0081 

.0138 

.0145 

10 

.2257 

.2308 

.2846 

80 

.0057 

.0101 

.0107 

15 

.0969 

.1129 

.1159 

90 

.0038 

.0074 

.0078 

20 

.0674 

.0821 

.0868 

100 

.0026 

.0054 

.0058 

4.2.2  Spherical  Immobile  Material  This  test  demonstrates  the  model’s  ability  to  simulate 
spherical  immobile  regions.  The  input  parameters  used  here  are  given  in  Table  4.4,  and  the  results 
of  model  executions  are  presented  in  Table  4.3. 

The  model  presented  here  agrees  favorably  with  the  analytic  solution  for  several  days,  then 
begins  to  show  a  markedly  higher  concentration.  By  the  end  of  the  test  period,  the  model  predicts 
about  twice  as  much  contaminant  at  the  well  as  the  analytic  solution.  This  is  due  to  the  way  the 
model  represents  spherical  immobile  material.  The  radius  of  a  typical  sphere,  and  the  number  of 

Table  4.4.  Input  Parameters  for  the  Spherical  Immobile  Material  Comparison  Test 


Qw 

0 

0m /e 

Dt 

H 

P 

f 

b 

.4, 

A'd 

1.16E-2 

.42 

2/3 

1.15E-10 

10. 

1.81E-3 

.4 

.05 

.5 

.0 
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divisions  in  the  immobile  region  is  input.  The  sphere  radius  is  divided  into  even  intervals,  with  the 
boundary  (the  mobile  region  concentration)  at  the  same  interval  from  the  boundary  of  the  spin  r<>. 
This  effectively  includes  the  extra  interval  in  the  immobile  element  and  since  the  extra  ‘shell'  adds 
a  substantial  volume  to  the  sphere,  the  amount  of  contaminant  released  by  desorption  is  much 
greater  than  expected.  For  the  simple  case  of  four  grid  points  in  a  sphere,  the  volume  increase  is 
95%!  As  the  number  of  grid  points  increases,  this  effect  is  diminished,  though  as  the  number  of 
grid  points  increases  in  the  immobile  region,  the  execution  time  increases  as  well. 

The  concentration  in  the  immobile  zone  decreases  more  rapidly  when  the  zone  is  spherical 
than  when  it  consists  of  layers.  This  is  because  for  the  same  volume,  the  spherical  aggregates  have 
a  much  larger  surface  area  in  contact  with  the  mobile  zone  than  the  layered  immobile  material. 
The  spherical  aggregates  have  six  times  as  much  surface  area  in  contact  with  the  boundary  than 
layered  immobile  material  of  the  same  volume  and  halfwidth.  For  instance,  with  one  cubic  meter 
and  a  halfwidth  of  0.05  meters,  the  surface  area  for  layered  immobile  material  is  10  square  meters, 
ami  for  spherical  immobile  material  is  60  square  meters.  This  also  means  that  the  contaminant 
concentration  in  the  spheres  decreases  more  rapidly  than  is  observed  for  layered  immobile  material 
when  all  else  is  identical. 

4-  J  Adveclion,  Dispersion,  and  Slow  Desorption  Tests 

The  comparison  undertaken  here  is  one  where  the  input  parameters  are  adjusted  so  that  the 
program  models  advection  only,  then  advection  and  dispersion,  and  finally  advection,  dispersion, 
and  slow  desorption.  The  temporal  progressions  for  these  three  model  executions  are  shown  in 
Figures  4.13,  4.14,  and  4.15,  with  the  concentrations  extracted  from  the  well  shown  in  Table  1.6. 
These  three  time  progressions  are  displayed  with  time  decreasing  with  the  perspective  (the  opposite 
of  other  similar  figures  in  this  chapter.)  The  input  parameters  are  shown  in  'Fable  4.5.  Advert  ion 
only  is  achieved  by  making  the  diffusion  in  the  immobile  region  so  small  it  becomes  negligible,  and 
dispersion  is  eliminated  by  setting  the  dispersivity  to  zero.  The  model  simulated  layered  immobile 
material  for  this  test  and  used  finite  elements  in  the  mobile  zone. 


Figure  4.13.  Mobile  Zone  Concentration  Profile,  Days  0-30,  Advection  Only  (reversed  perspec¬ 
tive) 

With  advection  only,  there  should  not  be  a  spreading  of  the  contaminant.  The  sharp  drop 
in  contaminant  should  be  advected  intact  toward  the  well  However,  because  of  the  representation 
in  the  model,  the  top  of  this  gradient  is  one  grid  point  closer  to  the  well  than  the  bottom  of  the 
gradient.  This,  coupled  with  the  slightly  increased  water  velocity  makes  the  top  of  this  gradient 
move  faster  than  the  bottom. 

The  model  equations  do  not  provide  smoothing  of  the  data  or  require  the  concentration  to  be 
positive  in  sign.  The  dispersion  induces  a  smoothing  of  the  data,  and  prevents  spurious  negative 
values  from  occurring.  Unfortunately,  when  the  model  simulates  advection  only,  the  shock  condition 
from  the  input  profile  causes  an  overshooting  that  yields  negative  concentrations.  When  the  input 
data  is  devoid  of  these  shock  conditions,  negative  concentrations  don't  appear. 

With  advection  only,  all  the  contaminant  would  be  removed  from  the  aquifer  When  the 
dispersion  coefficient  is  not  zero,  the  contaminant  is  spread  out  and  it  takes  longer  to  remove  the 
contaminant  from  the  aquifer.  By  comparing  the  final  profiles  in  Figures  4.14  a' .j  4.15  and  t he 
concentrations  extracted  from  the  well  (Table  4.6)  it  is  evident  that  the  slow  desorption  process  is 
especially  important  in  the  latter  part  of  the  remediation  effort. 
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Figure  4.14.  Mobile  Zone  Concentration  Profile,  Days  0-30,  Advection  and  Dispersion  (reversed 
perspective) 


Figure  4.15.  Mobile  Zone  Concentration  Profile,  Days  0-30,  Advection,  Dispersion,  and  Slow 
Desorption  (reversed  perspective) 


Table  4.5.  Input  Parameters  for  the  Advection,  Dispersion,  and  Slow  Desorption  Comparison 
Test 


Test 

Qw 

e 

dm/0 

D e 

H 

p 

f 

MM 

m 

warn 

Adv 

1.16E-2 

.42 

.5 

1.15E-20 

10. 

1.81E-3 

B 

mm 

Adv+Disp 

1.16E-2 

.42 

.5 

1.15E-20 

10. 

’  81E-3 

Kf 

npy 

All 

1.16  E-  2 

.42 

.5 

1 . 15E2- 10 

10. 

1.81E-3 

n 
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Table  4.6.  Arlvection,  Dispersion,  and  Slow  Desorption  Comparison  Test 


Day 

Adv 

Adv+Dis 

All 

Day 

Adv 

Adv+Dis 

All 

I  0 

1.000 

1.000 

1.000 

16 

.9454 

.7179 

.7845 

2 

.„C93 

1.000 

1.000 

IS 

.9268 

.6414 

.6520 

!  4 

1.000 

.9998 

.9998 

20 

.6582 

.5002 

.5164 

O 

.9995 

.9999 

.9999 

22 

.4670 

.3708 

.3929 

8 

.9988 

1.004 

1.004 

24 

.2854 

.2623 

.2900 

10 

1.004 

1.001 

1.001 

26 

.1357 

.1755 

.2084 

12 

1.013 

.9693 

.9698 

28 

.0401 

.1141 

.1511 

14 

1.004 

.8955 

.8980 

30 

-.001 

.0714 

.1117 

44  Pulse  Pumping 

The  first  test  above  had  a  duration  of  100  days.  In  this  section  looks  at  the  continuation  of 
that  test  out  to  400  days  with  different  pumping  rates.  In  this  first  test  (Pulse  Pumping  'lest  1) 
the  pumping  rate  was  one  cubic  meter  per  day  for  the  remaining  300  days  (see  Figure  4.16).  In  the 
next  test  (Puls-.  Pumping  Test  2)  the  pumping  rate  was  one  cubic  meter  per  day  for  days  100  to 
200.  then  at  1U00  cubic  meters  per  da\  (the  value  used  in  the  first  test  in  this  chapter)  from  day 
200  to  300,  and  then  back  to  one  cubic  meter  per  day  for  the  last  ICO  days  (see  Figure  4.17).  In 
the  final  example  (Pulse  Pumping  Test  3)  the  pumping  rate  was  50  cubic  meters  p  day  for  the 
entire  period  from  day  100  to  day  400  (see  Figure  4.18). 

As  can  be  seen  from  Figures  4.19  through  4.32,  if  the  pumping  is  ceased  before  the  level  in 
the  immobile  region  is  reduced  significantly,  the  concentrations  in  the  mobile  region  increase  with 
time.  This  is  not  an  unexpected  result,  and  has  already  been  observed  (22:633)  (19). 

The  creation  of  this  model  involved  one  assumption  about  diffusion:  dispersion  diffusion  in 
the  mobile  region,  and  therefore  diffusion  could  be  ignored.  With  no  pumping  there  is  no  guarantee 
that  diffusion  in  the  mobile  region  will  not  play  an  important  part  in  the  fate  of  the  pollutant  For 
this  reason,  a  small  pumping  rate  is  maintained,  even  if  it  represents  a  velocity  of  only  a  lew 
millimeters  per  day. 
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Figure  4.19.  Pulse  Pumping  Test  1;  Negligible  Flow,  Progress  from  Days  1UU-200 


Figure  4.21  Pulse  Pumping  Test  1,  Immobile  Region  Profile  at  Day  1100 


Figure  4.25.  Pulse  Pumping  Test  2;  Immobile  Region  Profile  at  Day  300 
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Figure  4.28.  Pulse  Pumping  Test  3;  Reduced  Flow,  Progress  from  Days  100-  200 
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\  .  Conclusions  and  Recommendations 

This  concluding  chapter  draws  together  the  research  presented  in  the  previous  chapters  It 
focuses  on  how  the  research  may  improve  the  ability  of  planners  to  estimate  the  needed  duration 
time  of  remediation  pumping.  It  also  shows  how  slow  sorption  and  desorption  of  chemicals  in 
porous  media  could  play  an  important  part  in  the  fate  of  pollutants.  Secondly,  it  enumerates  the 
conclusions  which  may  be  drawn  from  the  research.  Finally,  it  lists  recommendations  both  for 
improvements  and  for  the  follow-on  efforts  to  this  research. 

5.1  Summary 

The  literature  survey  revealed  that  research  in  the  last  few-  decades  has  contributed  to  an  im¬ 
proved  understanding  of  the  desorption  processes  in  contaminant  transfer  through  porous  media. 
Although  there’s  been  an  increase  in  our  abilities  to  mathematically  model  this  contaminant  trans¬ 
port,  there  exists  a  paucity  of  computer  based  solutions  that  incorporate  rate-limited  desorption. 

5.3  Conclusions 

Several  conclusions  could  be  drawn  from  the  comparison  of  this  computer  based  solut  ion  to  the 
contaminant  traits  oct.  equations  with  the  analytical  solution  provided  by  Goltz.  Unfortunately, 
without  a  large  database  for  comparison  with  field  data  or  laboratory  experiment,  we  can  not 
properly  evaluate  this  research.  An  initial  laboratory  experiment  will  be  conducted  shortly  by  Dr. 
Burris  (7),  but  not  soon  enough  to  add  to  this  thesis.  That  laboratory  experiment  will  consist  of 
monitoring  the  contaminant  levels  in  a  three  dimensional  physical  mode!  with  axial  symmetry.  A 
four  foot  diameter  tank  will  be  filled  with  a  sand  or  soil  material,  and  will  have  a  fully  penetrating 
well  in  the  center  of  the  tank.  The  contaminant  concentrations  at  the  well  and  in  the  tank  will  be 
monitored. 

If  the  laboratory  experiments  yield  favorable  comparisons  with  this  numerical  solution,  and 
display  similar  tailing  of  contaminants,  this  simple  model  could  be  distributed  to  planners  at  t>  >xw 
cleanup  sites  for  incorporation  into  their  computer  toolkit. 
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5.3  Recommendations 


Most  of  the  drawbacks  of  this  model  stem  from  the  assumptions  made  in  the  Assum/ilion .•> 
section  of  Chapter  1.  The  primary  emphasis  for  future  research  is  to  eliminate  these  assumptions 
(in  the  order  given): 


1.  Eliminate  the  necessity  for  spherical  symmetry;  make  the  mathematical  model  truly  two- 
dimensional. 

2.  Add  the  ability  to  simulate  multiple  wells,  including  injection  wells. 

3.  Create  a  user-friendly  input  port  for  the  contaminant  profile;  probably  based  on  mouse-driven 
inputs  on  an  Amiga  or  a  Sun  workstation. 

■1.  Add  the  ability  to  estimate  the  total  contaminant  remaining  in  the  aquifer,  and  the  rate  at 
which  it  will  desorb  at  normal  seepage  rates. 


Other  improvements  could  include  correcting  the  problem  with  the  apparent  radius  when 
using  diffusion  into  spherical  aggregates.  Only  experimentation  will  indicate  if  this  is  advised. 
Also,  instead  of  allowing  only  spherical  or  layered  immobile  material,  allow  the  third  possibility: 
cylindrical  immobile  material.  One  additional  enhancement  could  be  a  further  simplification  of  the 
existing  code  to  use  the  first-order  models  presented  by  Coats  ami  Smith  ( 10).  van  (lenuchten  (3C>). 
and  Valocchi  (33). 


5.J,  Remarks 

As  industrial  development  continues,  and  as  long  as  Murphy's  laws  apply,  there  will  be  unfor¬ 
tunate  instances  of  polluted  groundwater  to  address.  This,  together  with  the  need  to  cleanup  the 
existing  toxic  spill  sites,  make  this  an  important  piece  of  research  for  the  environmentally  conscious 


community.  This  research  has  endeavored  to  further  the  cause  of  numerical  modeling  applied  to 
porous  media,  and  to  provide  the  community  with  another  tool  to  help  in  the  field. 


Appendix  A.  Notation 


These  are  most  of  the  variables  or  constants  used  in  this  Thesis.  Parenthesized  names  are 
the  variable  names  used  in  MODEL1T  representing  the  same  data. 


Ai  —  [L]  Dispersivity  in  the  mobile  region  (alphaL) 

b  -  [L]  Characteristic  length  of  the  immobile  region  geometry,  the  halfwidth  or  radius  of  the 

immobile  region 

C  —  [A//L3]  Concentration  of  adsorbate  in  the  fluid  stream:  Solute  concentration 

Cm  —  [Af/L3]  Solute  concentration  vector  in  the  mobile  region  (Cm) 

C'i,n  —  [Af/L3]  Solute  concentration  vector  in  the  immobile  region  (Cim) 

C\m(l  —  [A// /^3]  Volume  averaged  solute  concentration  vector  for  the  immobile  region  (Cima) 

D  —  [ Z. 2 / T^]  Hydrodynamic  dispersion  coefficient 

Df  —  [L2/T\  Diffusion  coefficient  within  the  immobile  region  (De) 

D *  -  f L2/T]  Diffusion  coefficient  within  the  mobile  region 

/  [unitless]  Fraction  of  sorption  sites  in  direct  contact  with  mobile  water  (ef) 

//  —  [L]  Average  height  of  the  aquifer 

h  [T-1]  Exchange  rate  between  stagnant  and  mobile  water 
n’d  —  [L3/M]  Distribution  Coefficient  (Kd) 

Qu  —  [L3/T]  Flow  rate  at  the  well  (Flow),  extraction  being  positive 

r  ■  —  [L]  Radial  coordinate  within  the  mobile  region;  0  <  r  <  R 

R  —  [unitless]  Average  retardation  factor:  R  =  1  +  (retardation  is  a  friction-type  force  that 

0 

slows  the  progress  of  the  contaminant  relative  to  the  flow  of  groundwater  velocity,  typical 
values  arc  1 .  .  .33):  also  used  as  the  maximum  extent  of  the  contaminated  aquifer  in  the 
radial  dimension  [L] 

/?,„  -  [unitless]  Mobile  region  retardation  factor:  Rm  =  1  +  —  ^  -  —  (Rm) 


Hi 


[unitless]  Immobile  region  retardation  factor:  R„n  =  1  + 


(1  -f)  Ph’j 

Oim 


(Him) 


•S’  [A//L3]  Sorbed  solute  concentration,  or  [unitless]  sorbed  contaminant 

l  —  [T]  Time 

V  (r)  —  [Z/T]  Velocity  of  the  groundwater  through  the  interstices  of  the  aquifer  (mobile  regions 

only) 


-  Greek  — 


n 


.1 


[7~]  Basis  function  (temporal  extent) 

[ /.]  Spatial  coordinate  within  the  immobile  region 


[unit less]  Solute  capacity  ratio  of  the  immobile  to  mobile  regions:  0  —  - (beta) 


A-l 


0  --  [unitless]  Porosity  or  Total  water  content:  9  —  6m  +  #,m  (theta),  also  called  total  poro>ity 
(typical  values  for  granular  aquifers  are  20-40%  ) 

0,m  —  [unitless]  Immobile  water  content,  also  called  immobile  region  porosity 

Or,,  —  [unitless]  Mobile  water  content  (tlietaM),  also  called  mobile  region  porosity 

A  —  [unitless]  Fraction  representat  ion  of  stagnant  to  mobile  water 

n  —  [unitless]  Ratio  of  diameter  of  a  circle  to  its  circumference 

p  —  [M/L3]  Bulk  density  of  aquifer  material  (rho):  rho  is  also  used  as  a  constant  in  the  model 

r  —  [unitless]  Integration  variable 

0  —  [L]  Basis  functions  (spatial  extent),  also  [unitless]  ratio  of  mobile  to  total  water  content: 

4>  =  em/e  (phi) 
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Abstract 


This  paper  traces  the  development  of  the  understanding  of  contaminant  transport  in  an 
aquifer  for  a  radially  symmetric  region.  It  presents  the  progression  of  ideas  and  equations  lead¬ 
ing  to  the  equation  set  describing  the  advective/dispersive  mechanisms  coupled  with  rate-limited 
adsorption.  This  equation  set  is  converted  to  a  numerical  scheme  in  the  Fortran  language  using 
finite  elements  and  finite  differences.  The  resulting  model  is  tested  against  the  Laplace  transform 
solution  to  the  same  equation  set,  and  several  graphs  are  presented  detailing  the  comparison 
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